#define MECPlotting_cxx
#include "MECPlotting.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

#include <iostream>
#include <vector>
#include <fstream>
#include <sstream>
#include <string>
#include <list>
#include <functional>
#include <numeric>
#include <TClonesArray.h>
#include <TRefArray.h>
#include <TMath.h>

#include <TSystem.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TF1.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TTree.h>
#include "TDirectory.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TString.h"
#include <TLegend.h>
#include "TLine.h"


#define TRUE 0
#define FALSE 1
#define PROTON 2212
#define NEUTRON 2112
#define MUON 13


using namespace std;
using std::cout;
using std::endl;

TList* addHistogram1D(TList*, const char*, const char*, int, double, double);
TList* addHistogram2D(TList*, const char*, const char*, int, double, double, int, double, double);
TList* initialiseHists(TList*);
//int    fillHists(int, TChain*, TList*&);

//void MECPlotting::Loop()
int MECPlotting::analysis(int NEvents, TString filename, TString outSaveName)
{
  //Specify Error Flag
  int ISNOERROR = TRUE;  
  
  //Load input file
  ISNOERROR = MECPlotting::SelectFile(filename);
  if(ISNOERROR==FALSE)
  {
    cout << "Error: cannot read file: " << filename << endl;
    getchar();
  }
  
  //Initialise Histogram list
  TList* histList = new TList();

  //create hist templates  
  histList = initialiseHists(histList);

  ISNOERROR = fillHists(NEvents, histList);
  if(ISNOERROR==FALSE)
  {
    cout << "Error: cannot fill hists " << endl;
    getchar();
  }
  
  //Wrap Up:

  TFile histSave(outSaveName,"RECREATE");

  histList->Write();
  
  //Canvas creation and combined histograms goes here
    
  
  histSave.Close();


  if (ISNOERROR == FALSE)
  {
    cout << "Error: Exiting now ..." << endl;
    return(1);
  }
  else
  {
    cout << "Successful Run, Exiting now ..." << endl;
    return(0);
  }

}

int MECPlotting::fillHists(int NEvents, TList* &histList){
  
   //Masses are in MeV
   double M_n = 939.56563;
   double M_p = 938.272046;
   double M_m = 105.6583715;
   
         
   int reconMec = 0;
   int reconCcqe = 0;
   int reconRes = 0;
   int reconOther = 0;
   int recoMECBeforeCuts = 0;
   int recoCCQEBeforeCuts = 0;
   int recoRESBeforeCuts = 0;
   int recoOTHERBeforeCuts = 0;
   int recoALLBGBeforeCuts = 0;
   const int cutSize = 999;
   vector<int> mecAfterTmCut(cutSize, 0);
   vector<int> ccqeAfterTmCut(cutSize, 0);
   vector<int> resAfterTmCut(cutSize, 0);
   vector<int> otherAfterTmCut(cutSize, 0);
   vector<int> AllBGAfterTmCut(cutSize, 0);

   vector<double> deltaPhiPtCut(cutSize, 0);
   vector<double> targetMassCut(cutSize, 9000);
   TVector3 DeltaRecPt, DeltaSimPt;
   TVector3 RecKp_t, SimKp_t;
   TVector3 RecKmu_t, SimKmu_t;
   TVector3 nuDirRec, nuDirSim;

  
   if (fChain == 0) return 1;
   
   Long64_t nentries = 0;
   if(NEvents != 0) nentries = NEvents;
   else nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = MECPlotting::LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      std::cout << "Processing entry "<< jentry << std::endl;
      
      DeltaRecPt.SetXYZ(fDeltaRecPtfX,fDeltaRecPtfY,fDeltaRecPtfZ);      
      DeltaSimPt.SetXYZ(fDeltaSimPtfX,fDeltaSimPtfY,fDeltaSimPtfZ);  
      nuDirRec.SetXYZ(fNeutrinoRecfX,fNeutrinoRecfY,fNeutrinoRecfZ);
      nuDirRec = nuDirRec.Unit();
      nuDirSim.SetXYZ(fNeutrinoSimfPfX,fNeutrinoSimfPfY,fNeutrinoSimfPfZ);
      nuDirSim = nuDirSim.Unit();
      fProtonRecfP.SetXYZ(fProtonRecfPfX,fProtonRecfPfY,fProtonRecfPfZ);
      fProtonSimfP.SetXYZ(fProtonSimfPfX,fProtonSimfPfY,fProtonSimfPfZ);
      fMuonRecfP.SetXYZ(fMuonRecfPfX,fMuonRecfPfY,fMuonRecfPfZ);
      fMuonSimfP.SetXYZ(fMuonSimfPfX,fMuonSimfPfY,fMuonSimfPfZ);
      RecKp_t.SetXYZ(fRecKp_tfX, fRecKp_tfY, fRecKp_tfZ);
      SimKp_t.SetXYZ(fSimKp_tfX, fSimKp_tfY, fSimKp_tfZ);      
      RecKmu_t.SetXYZ(fRecKmu_tfX, fRecKmu_tfY, fRecKmu_tfZ);
      SimKmu_t.SetXYZ(fSimKmu_tfX, fSimKmu_tfY, fSimKmu_tfZ);
      
      int ISMEC = FALSE;
      int ISCCQE = FALSE;
      int ISCCRES = FALSE;
      int ISOTHER = FALSE;
      
      if(fNeutMode==2)
      {
        reconMec++;
        ISMEC = TRUE;
        std::cout << "Found Recon MEC , code is "<< fNeutMode << std::endl;
      }
      else if(fNeutMode==1)
      {
        reconCcqe++;
        ISCCQE = TRUE;
        std::cout << "Found Recon CCQE , code is " << fNeutMode << std::endl;
      }
      else if(fNeutMode==11  || fNeutMode==16)
      {
        reconRes++;
        ISCCRES = TRUE;
        std::cout << "Found Recon RES , code is " << fNeutMode << std::endl;
      }
      else{
        reconOther++;
        ISOTHER = TRUE;
      }
      
      if(ISMEC == TRUE){

        cout << "MEC nu Angle" << endl;
        ((TH1*)histList->FindObject("nuAngle"))->Fill(180*acos(nuDirRec.Z())/(TMath::Pi()));
        ((TH1*)histList->FindObject("nuAngle_t"))->Fill(180*acos(nuDirSim.Z())/(TMath::Pi()));
        
        //Proton Variables:       
         cout << "MEC Proton Momentum" << endl;
        ((TH1*)histList->FindObject("pPHistInFGD"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("pPHistInFGD_t"))->Fill(fProtonSimfP.Mag());

        cout << "MEC Proton Neutrino Angle Calc" << endl;
        double pnangle   = TMath::ACos(nuDirRec*fProtonRecfP/(nuDirRec.Mag()*fProtonRecfP.Mag()));
        double pnangle_t = TMath::ACos(nuDirSim*fProtonSimfP/(nuDirSim.Mag()*fProtonSimfP.Mag()));

         cout << "MEC Proton Neutrino Angle" << endl;
        ((TH1*)histList->FindObject("AngleHistInFGD"))->Fill(pnangle);
        ((TH1*)histList->FindObject("AngleHistInFGD_t"))->Fill(pnangle_t);

        cout << "MEC Proton Energy" << endl;
        ((TH1*)histList->FindObject("pEHistInFGD"))->Fill(fProtonRecfE);
        ((TH1*)histList->FindObject("pEHistInFGD_t"))->Fill(fProtonSimfE);
        
        //Muon Variables:
        cout << "MEC Muon Momentum" << endl;
        ((TH1*)histList->FindObject("mPHistInFGD"))->Fill(fMuonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistInFGD_t"))->Fill(fMuonSimfP.Mag());
       
        cout << "MEC Muon Neutrino Angle Calc" << endl;
        double mnangle =   TMath::ACos(nuDirRec*fMuonRecfP/(nuDirRec.Mag()*fMuonRecfP.Mag()));
        double mnangle_t = TMath::ACos(nuDirSim*fMuonSimfP/(nuDirSim.Mag()*fMuonSimfP.Mag()));
         
        cout << "MEC Muon Neutrino Angle" << endl;
        ((TH1*)histList->FindObject("mnAngleHistInFGD"))->Fill(mnangle);
        ((TH1*)histList->FindObject("mnAngleHistInFGD_t"))->Fill(mnangle_t); 
       
        cout << "MEC Muon Energy" << endl;        
        ((TH1*)histList->FindObject("mEHistInFGD"))->Fill(fMuonRecfE);
        ((TH1*)histList->FindObject("mEHistInFGD_t"))->Fill(fMuonSimfE);
        
        //Calculation Variables:
        
        ((TH1*)histList->FindObject("EnuRecHistFromP"))->Fill(fReconRecNuEnergy);
        ((TH1*)histList->FindObject("EnuRecHistFromP_t"))->Fill(fReconSimNuEnergy);
        ((TH1*)histList->FindObject("ReconTargetMassFromP_MEC"))->Fill(fReconRecTgtMass);
        ((TH1*)histList->FindObject("ReconTargetMassFromP_MEC_t"))->Fill(fReconSimTgtMass);

        ((TH1*)histList->FindObject("difP_tHistInFGD"))->Fill(DeltaRecPt.Mag());
        ((TH1*)histList->FindObject("difP_tHistInFGD_t"))->Fill(DeltaSimPt.Mag());

        ((TH1*)histList->FindObject("deltaPhiPt"))->Fill(fDeltaRecPhi_t);
        ((TH1*)histList->FindObject("deltaPhiPt_t"))->Fill(fDeltaSimPhi_t);
        
        ((TH1*)histList->FindObject("protonKt"))->Fill(RecKp_t.Mag());
        ((TH1*)histList->FindObject("protonKt_t"))->Fill(SimKp_t.Mag());
        
        ((TH1*)histList->FindObject("muonKt"))->Fill(RecKmu_t.Mag());
        ((TH1*)histList->FindObject("muonKt_t"))->Fill(SimKmu_t.Mag());
        
        ((TH2*)histList->FindObject("deltaPhiPtvsProtonKt"))->Fill(fDeltaRecPhi_t,RecKp_t.Mag());        
        ((TH2*)histList->FindObject("deltaPhiPtvsMuonKt"))->Fill(fDeltaRecPhi_t,RecKmu_t.Mag());
        ((TH2*)histList->FindObject("ProtonKtvsMuonKt"))->Fill(RecKp_t.Mag(),RecKmu_t.Mag());
        
        ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_MEC"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);        
        ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_MEC"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);  

        ((TH2*)histList->FindObject("EnuRecVsReconTargetMass_MEC"))->Fill(fReconRecNuEnergy,fReconRecTgtMass);
        ((TH2*)histList->FindObject("ProtonKtVsReconTargetMass_MEC"))->Fill(RecKp_t.Mag(),fReconRecTgtMass);
        ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_MEC"))->Fill(RecKmu_t.Mag(),fReconRecTgtMass);
        ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_MEC"))->Fill(RecKmu_t.Mag(),fReconRecTgtMass);

        
      }
      
      if(ISMEC == FALSE){
        
        cout << "Not MEC nu Angle " << 180*acos(nuDirRec.Z())/(TMath::Pi()) << endl;
        ((TH1*)histList->FindObject("nuAngle"))->Fill(180*acos(nuDirRec.Z())/(TMath::Pi()));
        ((TH1*)histList->FindObject("nuAngle_t"))->Fill(180*acos(nuDirSim.Z())/(TMath::Pi()));
        
        //Proton Variables:       
         cout << "Not MEC Proton Momentum " << fProtonRecfP.Mag() <<  endl;
        ((TH1*)histList->FindObject("pPHistInFGD_NotMEC"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("pPHistInFGD_NotMEC_t"))->Fill(fProtonSimfP.Mag());

        cout << "Not MEC Proton Neutrino Angle Calc" << endl;
        double pnangle   = TMath::ACos(nuDirRec*fProtonRecfP/(nuDirRec.Mag()*fProtonRecfP.Mag()));
        double pnangle_t = TMath::ACos(nuDirSim*fProtonSimfP/(nuDirSim.Mag()*fProtonSimfP.Mag()));

         cout << "Not MEC Proton Neutrino Angle " << pnangle << endl;
        ((TH1*)histList->FindObject("AngleHistInFGD_NotMEC"))->Fill(pnangle);
        ((TH1*)histList->FindObject("AngleHistInFGD_NotMEC_t"))->Fill(pnangle_t);

        cout << "Not MEC Proton Energy " << fProtonRecfE << endl;
        ((TH1*)histList->FindObject("pEHistInFGD_NotMEC"))->Fill(fProtonRecfE);
        ((TH1*)histList->FindObject("pEHistInFGD_NotMEC_t"))->Fill(fProtonSimfE);
        
        //Muon Variables:
        cout << "Not MEC Muon Momentum" << endl;
        ((TH1*)histList->FindObject("mPHistInFGD_NotMEC"))->Fill(fMuonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistInFGD_NotMEC_t"))->Fill(fMuonSimfP.Mag());
       
        cout << "Not MEC Muon Neutrino Angle Calc " << endl;
        double mnangle =   TMath::ACos(nuDirRec*fMuonRecfP/(nuDirRec.Mag()*fMuonRecfP.Mag()));
        double mnangle_t = TMath::ACos(nuDirSim*fMuonSimfP/(nuDirSim.Mag()*fMuonSimfP.Mag()));
         
        cout << "Not MEC Muon Neutrino Angle" << endl;
        ((TH1*)histList->FindObject("mnAngleHistInFGD_NotMEC"))->Fill(mnangle);
        ((TH1*)histList->FindObject("mnAngleHistInFGD_NotMEC_t"))->Fill(mnangle_t); 
       
        cout << "Not MEC Muon Energy " << fMuonRecfE << endl;
        ((TH1*)histList->FindObject("mEHistInFGD_NotMEC"))->Fill(fMuonRecfE);
        ((TH1*)histList->FindObject("mEHistInFGD_NotMEC_t"))->Fill(fMuonSimfE);
        
        //Calculation Variables:
         cout << "Not MEC Calc" << endl;     

        ((TH1*)histList->FindObject("EnuRecHistFromP_NotMEC"))->Fill(fReconRecNuEnergy);
        ((TH1*)histList->FindObject("EnuRecHistFromP_NotMEC_t"))->Fill(fReconSimNuEnergy);
        ((TH1*)histList->FindObject("ReconTargetMassFromP_NotMEC"))->Fill(fReconRecTgtMass);
        ((TH1*)histList->FindObject("ReconTargetMassFromP_NotMEC_t"))->Fill(fReconSimTgtMass);

        //The problem is here, fix it:
        ((TH1*)histList->FindObject("difP_tHistInFGD_NotMEC"))->Fill(DeltaRecPt.Mag());
        ((TH1*)histList->FindObject("difP_tHistInFGD_NotMEC_t"))->Fill(DeltaSimPt.Mag());

        ((TH1*)histList->FindObject("deltaPhiPt_NotMEC"))->Fill(fDeltaRecPhi_t);
        ((TH1*)histList->FindObject("deltaPhiPt_NotMEC_t"))->Fill(fDeltaSimPhi_t);
        
        ((TH1*)histList->FindObject("protonKt_NotMEC"))->Fill(RecKp_t.Mag());
        ((TH1*)histList->FindObject("protonKt_NotMEC_t"))->Fill(SimKp_t.Mag());
        
        ((TH1*)histList->FindObject("muonKt_NotMEC"))->Fill(RecKmu_t.Mag());
        ((TH1*)histList->FindObject("muonKt_NotMEC_t"))->Fill(SimKmu_t.Mag());
        
        
        ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_NotMEC"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);        
        ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_NotMEC"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);        
        
        ((TH2*)histList->FindObject("deltaPhiPtvsProtonKt_NotMEC"))->Fill(fDeltaRecPhi_t,RecKp_t.Mag());        
        ((TH2*)histList->FindObject("deltaPhiPtvsMuonKt_NotMEC"))->Fill(fDeltaRecPhi_t,RecKmu_t.Mag());
        ((TH2*)histList->FindObject("ProtonKtvsMuonKt_NotMEC"))->Fill(RecKp_t.Mag(),RecKmu_t.Mag());        

        ((TH2*)histList->FindObject("EnuRecVsReconTargetMass_NotMEC"))->Fill(fReconRecNuEnergy,fReconRecTgtMass); 
        ((TH2*)histList->FindObject("ProtonKtVsReconTargetMass_NotMEC"))->Fill(RecKp_t.Mag(),fReconRecTgtMass);
        ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_NotMEC"))->Fill(RecKmu_t.Mag(),fReconRecTgtMass);
        ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_NotMEC"))->Fill(RecKmu_t.Mag(),fReconRecTgtMass);
        
         cout << "Not MEC Finish" << endl;        

        
      }
      ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_Comb"))->Fill(RecKmu_t.Mag(),fReconRecTgtMass);    
      ((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_MECDivComb"))->Divide(((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_MEC")),((TH2*)histList->FindObject("MuonKtVsReconTargetMass_redBins_Comb")),1,1);
      ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_Comb"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);        
      ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_Comb"))->Fill(fDeltaRecPhi_t,fReconRecTgtMass);        
      ((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_MECDivComb"))->Divide(((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_MEC")),((TH2*)histList->FindObject("deltaPhiPtvsReconTargetMass_redBins_Comb")),1,1);        

      
      if(ISMEC == TRUE)
      {
        recoMECBeforeCuts++;
        ((TH1*)histList->FindObject("pPHistBeforeCuts_MEC"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistBeforeCuts_MEC"))->Fill(fMuonRecfP.Mag());
      }
      if(ISCCQE == TRUE)
      {
        recoCCQEBeforeCuts++;
        ((TH1*)histList->FindObject("pPHistBeforeCuts_CCQE"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistBeforeCuts_CCQE"))->Fill(fMuonRecfP.Mag());
        ((TH1*)histList->FindObject("ReconTargetMassFromP_CCQE"))->Fill(fReconRecTgtMass);
        ((TH1*)histList->FindObject("deltaPhiPt_CCQE"))->Fill(fDeltaRecPhi_t);
        ((TH1*)histList->FindObject("protonKt_CCQE"))->Fill(RecKp_t.Mag());
        ((TH1*)histList->FindObject("muonKt_CCQE"))->Fill(RecKmu_t.Mag());

      }     
      if(ISCCRES == TRUE)
      {
        recoRESBeforeCuts++;
        ((TH1*)histList->FindObject("pPHistBeforeCuts_RES"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistBeforeCuts_RES"))->Fill(fMuonRecfP.Mag());
        ((TH1*)histList->FindObject("ReconTargetMassFromP_RES"))->Fill(fReconRecTgtMass);
        ((TH1*)histList->FindObject("deltaPhiPt_RES"))->Fill(fDeltaRecPhi_t);
        ((TH1*)histList->FindObject("protonKt_RES"))->Fill(RecKp_t.Mag());
        ((TH1*)histList->FindObject("muonKt_RES"))->Fill(RecKmu_t.Mag());
      }
      if(ISOTHER == TRUE)
      {
        recoOTHERBeforeCuts++;
        ((TH1*)histList->FindObject("pPHistBeforeCuts_OTHER"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistBeforeCuts_OTHER"))->Fill(fMuonRecfP.Mag());
        ((TH1*)histList->FindObject("ReconTargetMassFromP_OTHER"))->Fill(fReconRecTgtMass);
        ((TH1*)histList->FindObject("deltaPhiPt_OTHER"))->Fill(fDeltaRecPhi_t);
        ((TH1*)histList->FindObject("protonKt_OTHER"))->Fill(RecKp_t.Mag());
        ((TH1*)histList->FindObject("muonKt_OTHER"))->Fill(RecKmu_t.Mag());
      }
      if(ISMEC == FALSE)
      {
        recoALLBGBeforeCuts++;
        ((TH1*)histList->FindObject("pPHistBeforeCuts_NotMEC"))->Fill(fProtonRecfP.Mag());
        ((TH1*)histList->FindObject("mPHistBeforeCuts_NotMEC"))->Fill(fMuonRecfP.Mag());
      }
      
      //Cuts for feesability:
      cout << "Cut loop" << endl;        

      
      double phiPtCutNow = 0.0;
      double tgtMassCutNow = 0.6;
      for(double i = 0; i < 20; i++)
      {
        for(double j = 0; j < 20; j++)
        {
          //int ele = j + 20*i;
          deltaPhiPtCut[j+(20*i)] = phiPtCutNow;
          targetMassCut[j+(20*i)] = tgtMassCutNow;
          tgtMassCutNow = tgtMassCutNow + 0.05;
        }
        tgtMassCutNow = 0.6;
        phiPtCutNow = phiPtCutNow + 0.1;
      }

      deltaPhiPtCut[0] = 0.4;
      targetMassCut[0] = 0.95;
      
      for(int cutNum = 0; cutNum<cutSize; cutNum++)
      {
        if(fDeltaRecPhi_t > deltaPhiPtCut[cutNum] && fReconRecTgtMass < targetMassCut[cutNum])
        {
          if(ISMEC == TRUE) mecAfterTmCut[cutNum]++; 
          if(ISCCQE == TRUE) ccqeAfterTmCut[cutNum]++;
          if(ISCCRES == TRUE) resAfterTmCut[cutNum]++; 
          if(ISOTHER == TRUE) otherAfterTmCut[cutNum]++; 
          if(ISMEC == FALSE) AllBGAfterTmCut[cutNum]++; 
        }
      }
           
      cout << "Apply Cuts " << endl;     
      if(fDeltaRecPhi_t > deltaPhiPtCut[0] && fReconRecTgtMass < targetMassCut[0])
      {
        if(ISMEC == TRUE)
        {
          ((TH1*)histList->FindObject("pPHistAfterCuts_MEC"))->Fill(fProtonRecfP.Mag());
          ((TH1*)histList->FindObject("mPHistAfterCuts_MEC"))->Fill(fMuonRecfP.Mag());
        }
        if(ISCCQE == TRUE)
        {
          ((TH1*)histList->FindObject("pPHistAfterCuts_CCQE"))->Fill(fProtonRecfP.Mag());
          ((TH1*)histList->FindObject("mPHistAfterCuts_CCQE"))->Fill(fMuonRecfP.Mag());
        }     
        if(ISCCRES == TRUE)
        {
          ((TH1*)histList->FindObject("pPHistAfterCuts_RES"))->Fill(fProtonRecfP.Mag());
          ((TH1*)histList->FindObject("mPHistAfterCuts_RES"))->Fill(fMuonRecfP.Mag());
        }
        if(ISOTHER == TRUE)
        {
          ((TH1*)histList->FindObject("pPHistAfterCuts_OTHER"))->Fill(fProtonRecfP.Mag());
          ((TH1*)histList->FindObject("mPHistAfterCuts_OTHER"))->Fill(fMuonRecfP.Mag());
        }
        if(ISMEC == FALSE)
        {
          ((TH1*)histList->FindObject("pPHistAfterCuts_NotMEC"))->Fill(fProtonRecfP.Mag());
          ((TH1*)histList->FindObject("mPHistAfterCuts_NotMEC"))->Fill(fMuonRecfP.Mag());
        }
      }

   } // End of Event loop
   
  vector<double> pur(cutSize, 0);
  vector<double> eff(cutSize, 0);
  vector<double> eff2(cutSize, 0);
  
  
  //double purA, effA, purB, effB, purC, effC, purD, effD;
  
  cout << "Variables after Tm cut: " << endl;
  cout << " Staistic          || deltaPhiPtCut || targetMassCut || Value" << endl;
  cout << "-------------------||---------------||---------------||---------" << endl;
  for(int ve = 0; ve <cutSize; ve++){
    cout.precision(2);
    cout << std::fixed;
    if(ve < 10) cout << " mecAfterTmCut" << ve << "    || " << deltaPhiPtCut[ve] << "          || " << targetMassCut[ve] << "       || "<<  mecAfterTmCut[ve] << endl;
    else        cout << " mecAfterTmCut" << ve <<  "   || " << deltaPhiPtCut[ve] << "          || " << targetMassCut[ve] << "       || "<<  mecAfterTmCut[ve] << endl;
    if(deltaPhiPtCut[ve] == 0 && targetMassCut[ve] == 9000 ) break;
  }
  cout << "-------------------||---------------||---------------||---------" << endl;
  cout << " reconMec          || " << reconMec << endl;
  cout << " reconCcqe         || " << reconCcqe << endl << endl << endl;
  
  for(int c = 0; c <cutSize; c++){
    pur[c] = 100*(double)mecAfterTmCut[c]/(mecAfterTmCut[c]+AllBGAfterTmCut[c]);
    eff[c] = 100*(double)mecAfterTmCut[c]/reconMec;
    eff2[c] = 100*(double)mecAfterTmCut[c]/recoMECBeforeCuts;

  } 
  
  cout << "Statistics after Tm cut: " << endl;
  cout << " CutNum || deltaPhiPtCut || targetMassCut || Purity  || Effeciency || Effeciency2" << endl;
  cout << "--------||---------------||---------------||---------||------------|| -------------- " << endl;
  for(int s = 0; s <cutSize; s++){
    cout.precision(2);
    cout << std::fixed;
    if(s < 10) cout << "  " << s << "     ||  " << deltaPhiPtCut[s] << "         || " << targetMassCut[s] << "       || " << pur[s]  << "   || " << eff[s] << "    || " << eff2[s]<< endl;
    else       cout << "  " << s <<  "    ||  " << deltaPhiPtCut[s] << "         || " << targetMassCut[s] << "       || " << pur[s]  << "   || " << eff[s] << "    || " << eff2[s]<< endl;
    if(deltaPhiPtCut[s] == 0 && targetMassCut[s] == 9000 ) break;
  }
  cout << "--------||---------------||---------------||---------||--------------- " << endl << endl << endl;
  
  ofstream statsOut;
  statsOut.open("statsOut.txt");
  
  statsOut << "Variables after Tm cut: " << endl;
  statsOut << " Staistic          || deltaPhiPtCut || targetMassCut || Value" << endl;
  statsOut << "-------------------||---------------||---------------||---------" << endl;
  for(int ve = 0; ve <cutSize; ve++){
    statsOut.precision(2);
    statsOut << std::fixed;
    if(ve < 10) statsOut << " mecAfterTmCut" << ve << "    || " << deltaPhiPtCut[ve] << "          || " << targetMassCut[ve] << "       || "<<  mecAfterTmCut[ve] << endl;
    else        statsOut << " mecAfterTmCut" << ve <<  "   || " << deltaPhiPtCut[ve] << "          || " << targetMassCut[ve] << "       || "<<  mecAfterTmCut[ve] << endl;
    if(deltaPhiPtCut[ve] == 0 && targetMassCut[ve] == 9000 ) break;
  }
  statsOut << "-------------------||---------------||---------------||---------" << endl;
  statsOut << " reconMec          || " << reconMec << endl;
  statsOut << " reconCcqe         || " << reconCcqe << endl << endl << endl;
  
  statsOut << "Statistics after Tm cut: " << endl;
  statsOut << " CutNum || deltaPhiPtCut || targetMassCut || Purity  ||| Effeciency || Effeciency2 " << endl;
  statsOut << "--------||---------------||---------------||---------||-------------||---------------- " << endl;
  for(int s = 0; s <cutSize; s++){
    statsOut.precision(2);
    statsOut << std::fixed;
    if(s < 10) statsOut << "  " << s << "     ||  " << deltaPhiPtCut[s] << "         || " << targetMassCut[s] << "       || " << pur[s]  << "   || " << eff[s] << "    || " << eff2[s]<< endl;
    else       statsOut << "  " << s <<  "    ||  " << deltaPhiPtCut[s] << "         || " << targetMassCut[s] << "       || " << pur[s]  << "   || " << eff[s] << "    || " << eff2[s]<< endl;
    if(deltaPhiPtCut[s] == 0 && targetMassCut[s] == 9000 ) break;
  }
  statsOut << "--------||---------------||---------------||---------||--------------- " << endl << endl << endl;

  statsOut.close();
  
  ofstream statsOutCSV;
  statsOutCSV.open("statsOut.csv");
  
  for(int s = 0; s <cutSize; s++){
    statsOutCSV.precision(2);
    statsOutCSV << std::fixed;
    statsOutCSV << s << ", " << deltaPhiPtCut[s] << ", " << targetMassCut[s] << ", " << pur[s] << ", " << eff[s] << ", " << eff2[s] <<  endl;
  }
  
  statsOutCSV.close();
  
  return 0;  
}



TList* addHistogram1D(TList* histList, const char* name, const char* title, int nbins, double min,double max){
  TH1* hist = new TH1D(name, title, nbins, min, max);
  histList->Add(hist);
  return(histList);
}



TList* addHistogram2D(TList* histList, const char* name, const char* title, int nbins_x, double min_x, double max_x, int nbins_y, double min_y, double max_y){
  TH2* hist = new TH2D(name, title, nbins_x, min_x, max_x, nbins_y, min_y, max_y);
  histList->Add(hist);
  return(histList);
}

TList* initialiseHists(TList* histList){
  
  //histList = addHistogram1D(histList, "FGDnuEHist","Histogram of Neutrino Energy inside FGD",300,0,14000);
  histList = addHistogram1D(histList, "pEHistInFGD","Histogram of MEC Proton Energy (final final)",50,0.8,2);
  histList = addHistogram1D(histList, "pPHistInFGD","Histogram of MEC Proton Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "AngleHistInFGD","Histogram of MEC Proton Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "mEHistInFGD","Histogram of MEC Muon Energy (final final)",100,0,2);
  histList = addHistogram1D(histList, "mPHistInFGD","Histogram of MEC Muon Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mnAngleHistInFGD","Histogram of MEC Muon to Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "pEHistInFGD_NotMEC","Histogram of Not MEC Proton Energy (final final)",50,0.8,2);
  histList = addHistogram1D(histList, "pPHistInFGD_NotMEC","Histogram of Not MEC Proton Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "AngleHistInFGD_NotMEC","Histogram of MEC Proton Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "mEHistInFGD_NotMEC","Histogram of Not MEC Muon Energy (final final)",100,0,2);
  histList = addHistogram1D(histList, "mPHistInFGD_NotMEC","Histogram of Not MEC Muon Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mnAngleHistInFGD_NotMEC","Histogram of Not MEC Muon to Neutrino angle (final final)",50,0,3);
  
  histList = addHistogram1D(histList, "pEHistInFGD_t","Histogram of MEC Proton Energy (final final)",50,0.8,2);
  histList = addHistogram1D(histList, "pPHistInFGD_t","Histogram of MEC Proton Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "AngleHistInFGD_t","Histogram of MEC Proton Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "mEHistInFGD_t","Histogram of MEC Muon Energy (final final)",100,0,2);
  histList = addHistogram1D(histList, "mPHistInFGD_t","Histogram of MEC Muon Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mnAngleHistInFGD_t","Histogram of MEC Muon to Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "pEHistInFGD_NotMEC_t","Histogram of Not MEC Proton Energy (final final)",50,8,2);
  histList = addHistogram1D(histList, "pPHistInFGD_NotMEC_t","Histogram of Not MEC Proton Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "AngleHistInFGD_NotMEC_t","Histogram of MEC Proton Neutrino angle (final final)",50,0,3);
  histList = addHistogram1D(histList, "mEHistInFGD_NotMEC_t","Histogram of Not MEC Muon Energy (final final)",100,0,2);
  histList = addHistogram1D(histList, "mPHistInFGD_NotMEC_t","Histogram of Not MEC Muon Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mnAngleHistInFGD_NotMEC_t","Histogram of Not MEC Muon to Neutrino angle (final final)",50,0,3);
  
  histList = addHistogram1D(histList, "EnuRecHist","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14);
  histList = addHistogram1D(histList, "EnuRecHist_t","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14);
  histList = addHistogram1D(histList, "ReconTargetMass_MEC","Histogram of Recontructed MEC Target Mass inside FGD",600,-2,6);
  histList = addHistogram1D(histList, "ReconTargetMass_MEC_t","Histogram of Recontructed MEC Target Mass inside FGD",600,-2,6);
  histList = addHistogram1D(histList, "EnuRecHist_NotMEC","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14);
  histList = addHistogram1D(histList, "EnuRecHist_NotMEC_t","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14);
  histList = addHistogram1D(histList, "ReconTargetMass_NotMEC","Histogram of Recontructed not MEC Target Mass inside FGD",600,-2,6);
  histList = addHistogram1D(histList, "ReconTargetMass_NotMEC_t","Histogram of Recontructed not MEC Target Mass inside FGD",600,-2,6);

  
  histList = addHistogram1D(histList, "EnuRecHistFromP","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14);
  histList = addHistogram1D(histList, "EnuRecHistFromP_t","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_MEC_t","Histogram of Recontructed MEC Target Mass inside FGD",600,-2,6);
  histList = addHistogram1D(histList, "EnuRecHistFromP_NotMEC","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14);
  histList = addHistogram1D(histList, "EnuRecHistFromP_NotMEC_t","Histogram of Recontructed Neutrino Energy inside FGD (Not MEC) ",300,0,14);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_NotMEC_t","Histogram of Recontructed not MEC Target Mass inside FGD",600,-2,6);

  histList = addHistogram1D(histList, "ReconTargetMassFromP_MEC","Histogram of Recontructed MEC Target Mass inside FGD",300,-0.5,3);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_CCQE","Histogram of Recontructed CCQE Target Mass inside FGD",300,-0.5,3);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_RES","Histogram of Recontructed RES Target Mass inside FGD",300,-0.5,3);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_OTHER","Histogram of Recontructed OTHER Target Mass inside FGD",300,-0.5,3);
  histList = addHistogram1D(histList, "ReconTargetMassFromP_NotMEC","Histogram of Recontructed not MEC Target Mass inside FGD",300,-0.5,3);



  histList = addHistogram1D(histList, "EnuRecHist_comb","Histogram of Recontructed Neutrino Energy inside FGD (MEC) ",300,0,14);
  histList = addHistogram1D(histList, "ReconTargetMass_comb","Histogram of Recontructed MEC Target Mass inside FGD",600,-2,6);


  histList = addHistogram1D(histList, "nOfFinalStateNucleons_MEC","Number of Nucleons in final state of MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfFinalStateProtons_MEC","Number of Protons in final state of MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfFinalStateNeutrons_MEC","Number of Neutrons in final state of MEC event",10,0,10);

  histList = addHistogram1D(histList, "nPIDS_MEC","Number of recon PIDS in MEC Event",10,0,10);
  histList = addHistogram1D(histList, "nOfReconProtons_MEC","Number of recon Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfReconMuons_MEC","Number of recon Muons in MEC event",10,0,10);

  histList = addHistogram1D(histList, "nOfTruthProtons_MEC","Number of true Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfTruthMuons_MEC","Number of true Muons in MEC event",10,0,10);

  histList = addHistogram1D(histList, "nOfTruthProtonsTPC_MEC","Number of true Protons in MEC event",10,0,10);
  histList = addHistogram1D(histList, "nOfTruthMuonsTPC_MEC","Number of true Muons in MEC event",10,0,10);

  histList = addHistogram2D(histList, "nOfTrueProtonsVsReconProtons","nOfTrueProtonsVsReconProtons",10,0,10,10,0,10);
  histList = addHistogram2D(histList, "nOfTrueProtonsTPCVsReconProtons","nOfTrueProtonsVsReconProtons",10,0,10,10,0,10);
  histList = addHistogram2D(histList, "nOfTrueProtonsVsProtonMom","nOfTrueProtonsVsReconProtonMom",10,0,10,50,0,12);

  histList = addHistogram2D(histList, "EnuRecVsReconTargetMass_MEC","Histogram of Recontructed Neutrino Energy vs Recon Target Mass (MEC) ",300,0,14,600,-4,4);
  histList = addHistogram2D(histList, "EnuRecVsReconTargetMass_NotMEC","Histogram of Recontructed Neutrino Energy vs Recon Target Mass (Not MEC) ",300,0,14,600,-4,4);

  histList = addHistogram1D(histList, "nuAngle","Histogram of incoming neutrino angle",100,0,3);
  histList = addHistogram1D(histList, "nuAngle_t","Histogram of incoming neutrino angle",100,0,3);


  histList = addHistogram1D(histList, "pP_tHistInFGD","Histogram of MEC Proton Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "pP_tHistInFGD_NotMEC","Histogram of Not MEC Proton Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mP_tHistInFGD","Histogram of MEC muon Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mP_tHistInFGD_NotMEC","Histogram of Not MEC muon Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "difP_tHistInFGD","Histogram of MEC dif between muon and proton Pt",50,-2,2);
  histList = addHistogram1D(histList, "difP_tHistInFGD_NotMEC","Histogram of Not MEC dif between muon and proton Pt",50,-2,2);

  histList = addHistogram1D(histList, "pP_tHistInFGD_t","Histogram of MEC Proton Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "pP_tHistInFGD_NotMEC_t","Histogram of Not MEC Proton Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mP_tHistInFGD_t","Histogram of MEC muon Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "mP_tHistInFGD_NotMEC_t","Histogram of Not MEC muon Transverse Momentum (final final)",50,-1,3);
  histList = addHistogram1D(histList, "difP_tHistInFGD_t","Histogram of MEC dif between muon and proton Pt",50,-2,2);
  histList = addHistogram1D(histList, "difP_tHistInFGD_NotMEC_t","Histogram of Not MEC dif between muon and proton Pt",50,-2,2);

  histList = addHistogram1D(histList, "deltaPhiPt","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_CCQE","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_RES","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_OTHER","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_NotMEC","Histogram of deltaPhiPt",100,0,3);
  
  histList = addHistogram1D(histList, "deltaPhiPt_t","Histogram of deltaPhiPt",100,0,3);
  histList = addHistogram1D(histList, "deltaPhiPt_NotMEC_t","Histogram of deltaPhiPt",100,0,3);
  
  histList = addHistogram1D(histList, "protonKt","Histogram of MEC protonKt",150,0,1.1);
  histList = addHistogram1D(histList, "protonKt_CCQE","Histogram of CCCQE protonKt",150,0,1.1);
  histList = addHistogram1D(histList, "protonKt_RES","Histogram of RES protonKt",150,0,1.1);
  histList = addHistogram1D(histList, "protonKt_OTHER","Histogram of OTHER protonKt",150,0,1.1);
  histList = addHistogram1D(histList, "protonKt_NotMEC","Histogram of Not MEC protonKt",150,0,1.1);
  
  histList = addHistogram2D(histList, "deltaPhiPtvsProtonKt","Histogram of deltaPhiPtvsProtonKt",100,0,3,100,0,2);
  histList = addHistogram2D(histList, "deltaPhiPtvsMuonKt","Histogram of deltaPhiPtvsProtonKt",100,0,3,100,0,2);
  histList = addHistogram2D(histList, "ProtonKtvsMuonKt","Histogram of deltaPhiPtvsProtonKt",100,0,2,100,0,2); 

  histList = addHistogram2D(histList, "deltaPhiPtvsProtonKt_NotMEC","Histogram of deltaPhiPtvsProtonKt",100,0,3,100,0,2);
  histList = addHistogram2D(histList, "deltaPhiPtvsMuonKt_NotMEC","Histogram of deltaPhiPtvsProtonKt",100,0,3,100,0,2);
  histList = addHistogram2D(histList, "ProtonKtvsMuonKt_NotMEC","Histogram of deltaPhiPtvsProtonKt",100,0,2,100,0,2); 
  
  histList = addHistogram2D(histList, "ProtonKtVsReconTargetMass_NotMEC","Histogram of ProtonKtVsReconTargetMass_NotMEC",300,0,2,300,-4,4); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_NotMEC","Histogram of MuonKtVsReconTargetMass_NotMEC",300,0,2,300,-4,4); 
  histList = addHistogram2D(histList, "ProtonKtVsReconTargetMass_MEC","Histogram of ProtonKtVsReconTargetMass_MEC",300,0,2,300,-4,4); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_MEC","Histogram of MuonKtVsReconTargetMass_MEC",300,0,2,300,-4,4); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_redBins_MEC","Histogram of MuonKtVsReconTargetMass_MEC RedBins",10,0,0.55,10,0.1,1.3); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_redBins_NotMEC","Histogram of MuonKtVsReconTargetMass_NotMEC RedBins",10,0,0.55,10,0.1,1.3); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_redBins_Comb","Histogram of MuonKtVsReconTargetMass_NotMEC RedBins",10,0,0.55,10,0.1,1.3); 
  histList = addHistogram2D(histList, "MuonKtVsReconTargetMass_redBins_MECDivComb","Histogram of MuonKtVsReconTargetMass MEC contribution",10,0,0.55,10,0.1,1.3); 

  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_MEC","Histogram of deltaPhiPtvsReconTargetMass_MEC",100,0,3,300,-4,4); 
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_NotMEC","Histogram of deltaPhiPtvsReconTargetMass_NotMEC",100,0,3,300,-4,4); 
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_Comb","Histogram of deltaPhiPtvsReconTargetMass_Comb",100,0,3,300,-4,4); 
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_redBins_Comb","Histogram of deltaPhiPtvsReconTargetMass_Comb",25,0,3,20,0,1.4);
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_redBins_MEC","Histogram of deltaPhiPtvsReconTargetMass_redBins_MEC",25,0,3,20,0,1.4);
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_redBins_NotMEC","Histogram of deltaPhiPtvsReconTargetMass_redBins_NotMEC",25,0,3,20,0,1.4);
  histList = addHistogram2D(histList, "deltaPhiPtvsReconTargetMass_redBins_MECDivComb","Histogram of deltaPhiPtvsReconTargetMass_redBins_MECDivComb",25,0,3,20,0,1.4);
  
  histList = addHistogram1D(histList, "protonKt_t","Histogram of MEC protonKt",150,0,1.1);
  histList = addHistogram1D(histList, "protonKt_NotMEC_t","Histogram of Not MEC protonKt",150,0,1.1);

  histList = addHistogram1D(histList, "muonKt","Histogram of MEC muonKt",150,0,1.1);
  histList = addHistogram1D(histList, "muonKt_CCQE","Histogram of CCCQE muonKt",150,0,1.1);
  histList = addHistogram1D(histList, "muonKt_RES","Histogram of RES muonKt",150,0,1.1);
  histList = addHistogram1D(histList, "muonKt_OTHER","Histogram of OTHER muonKt",150,0,1.1);
  histList = addHistogram1D(histList, "muonKt_NotMEC","Histogram of Not MEC muonKt",150,0,1.1);
  
  histList = addHistogram1D(histList, "muonKt_t","Histogram of MEC muonKt",150,0,1.1);
  histList = addHistogram1D(histList, "muonKt_NotMEC_t","Histogram of Not MEC muonKt",150,0,1.1);  


  histList = addHistogram1D(histList, "pPHistAfterCuts_MEC","Histogram of  MEC Proton Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistAfterCuts_CCQE","Histogram of  CCQE Proton Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistAfterCuts_RES","Histogram of  CCRES Proton Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistAfterCuts_OTHER","Histogram of  OTHER Proton Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistAfterCuts_NotMEC","Histogram of  Not MEC Proton Momentum after cuts",70,0,3);

  histList = addHistogram1D(histList, "pPHistBeforeCuts_MEC","Histogram of  MEC Proton Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_CCQE","Histogram of  CCQE Proton Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_RES","Histogram of  CCRES Proton Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_OTHER","Histogram of  OTHER Proton Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "pPHistBeforeCuts_NotMEC","Histogram of  Not MEC Proton Momentum before cuts",70,0,3);

  histList = addHistogram1D(histList, "mPHistAfterCuts_MEC","Histogram of  MEC Muon Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistAfterCuts_CCQE","Histogram of  CCQE Muon Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistAfterCuts_RES","Histogram of  CCRES Muon Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistAfterCuts_OTHER","Histogram of  OTHER Muon Momentum after cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistAfterCuts_NotMEC","Histogram of  Not MEC Muon Momentum after cuts",70,0,3);

  histList = addHistogram1D(histList, "mPHistBeforeCuts_MEC","Histogram of  MEC Muon Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_CCQE","Histogram of  CCQE Muon Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_RES","Histogram of  CCRES Muon Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_OTHER","Histogram of  OTHER Muon Momentum before cuts",70,0,3);
  histList = addHistogram1D(histList, "mPHistBeforeCuts_NotMEC","Histogram of  Not MEC Muon Momentum before cuts",70,0,3);

  return(histList);
 
}



/*
  // Declaration of leaf types
 TLorentzVector  *fMuonRec;
   TVector3        fMuonRecfP;
   Double_t        fMuonRecfE;
 TVector3        *fMuonRecPt;
 //TVector3        fMuonRecPt;
   Double_t        fMuonRecPtfX;
   Double_t        fMuonRecPtfY;
   Double_t        fMuonRecPtfZ;
   Double_t        fMuonRecPtMag;
   Double_t        fMuonRecAlpha;
   Double_t        fMuonRecQ2;
 TVector3        *fMuonRecVertex;
   Double_t        fMuonRecVertexfX;
   Double_t        fMuonRecVertexfY;
   Double_t        fMuonRecVertexfZ;
 TVector3        *fMuonRecEndPos;
   Double_t        fMuonRecEndPosfX;
   Double_t        fMuonRecEndPosfY;
   Double_t        fMuonRecEndPosfZ;
 TVector3        *fMuonRecFlightPath;
   Double_t        fMuonRecFlightPathfX;
   Double_t        fMuonRecFlightPathfY;
   Double_t        fMuonRecFlightPathfZ;
   Double_t        fMuonRecChi2;
   Int_t           fMuonRecNDOF;
   Int_t           fMuonRecNdets[2];
   Int_t           fMuonRecNhits[3];
   Int_t           fMuonRecCharge[3];
   Int_t           fMuonRecID;
   Int_t           fMuonRecPDG;
   Int_t           fMuonRecType;
 TLorentzVector  *fMuonSim;
   TVector3        fMuonSimfP;
   Double_t        fMuonSimfE;
 TVector3        *fMuonSimPt;
   Double_t        fMuonSimPtfX;
   Double_t        fMuonSimPtfY;
   Double_t        fMuonSimPtfZ;
   Double_t        fMuonSimPtMag;
   Double_t        fMuonSimAlpha;
   Double_t        fMuonSimQ2;
 TVector3        *fMuonSimVertex;
   Double_t        fMuonSimVertexfX;
   Double_t        fMuonSimVertexfY;
   Double_t        fMuonSimVertexfZ;
 TVector3        *fMuonSimEndPos;
   Double_t        fMuonSimEndPosfX;
   Double_t        fMuonSimEndPosfY;
   Double_t        fMuonSimEndPosfZ;
 TVector3        *fMuonSimFlightPath;
   Double_t        fMuonSimFlightPathfX;
   Double_t        fMuonSimFlightPathfY;
   Double_t        fMuonSimFlightPathfZ;
   Double_t        fMuonSimChi2;
   Int_t           fMuonSimNDOF;
   Int_t           fMuonSimNdets[2];
   Int_t           fMuonSimNhits[3];
   Int_t           fMuonSimCharge[3];
   Int_t           fMuonSimID;
   Int_t           fMuonSimPDG;
   Int_t           fMuonSimType;
 TLorentzVector  *fProtonRec;
   TVector3        fProtonRecfP;
   Double_t        fProtonRecfE;
 TVector3        *fProtonRecPt;
   Double_t        fProtonRecPtfX;
   Double_t        fProtonRecPtfY;
   Double_t        fProtonRecPtfZ;
   Double_t        fProtonRecPtMag;
   Double_t        fProtonRecAlpha;
   Double_t        fProtonRecQ2;
 TVector3        *fProtonRecVertex;
   Double_t        fProtonRecVertexfX;
   Double_t        fProtonRecVertexfY;
   Double_t        fProtonRecVertexfZ;
 TVector3        *fProtonRecEndPos;
   Double_t        fProtonRecEndPosfX;
   Double_t        fProtonRecEndPosfY;
   Double_t        fProtonRecEndPosfZ;
 TVector3        *fProtonRecFlightPath;
   Double_t        fProtonRecFlightPathfX;
   Double_t        fProtonRecFlightPathfY;
   Double_t        fProtonRecFlightPathfZ;
   Double_t        fProtonRecChi2;
   Int_t           fProtonRecNDOF;
   Int_t           fProtonRecNdets[2];
   Int_t           fProtonRecNhits[3];
   Int_t           fProtonRecCharge[3];
   Int_t           fProtonRecID;
   Int_t           fProtonRecPDG;
   Int_t           fProtonRecType;
 TLorentzVector  *fProtonSim;
   TVector3        fProtonSimfP;
   Double_t        fProtonSimfE;
 TVector3        *fProtonSimPt;
   Double_t        fProtonSimPtfX;
   Double_t        fProtonSimPtfY;
   Double_t        fProtonSimPtfZ;
   Double_t        fProtonSimPtMag;
   Double_t        fProtonSimAlpha;
   Double_t        fProtonSimQ2;
 TVector3        *fProtonSimVertex;
   Double_t        fProtonSimVertexfX;
   Double_t        fProtonSimVertexfY;
   Double_t        fProtonSimVertexfZ;
 TVector3        *fProtonSimEndPos;
   Double_t        fProtonSimEndPosfX;
   Double_t        fProtonSimEndPosfY;
   Double_t        fProtonSimEndPosfZ;
 TVector3        *fProtonSimFlightPath;
   Double_t        fProtonSimFlightPathfX;
   Double_t        fProtonSimFlightPathfY;
   Double_t        fProtonSimFlightPathfZ;
   Double_t        fProtonSimChi2;
   Int_t           fProtonSimNDOF;
   Int_t           fProtonSimNdets[2];
   Int_t           fProtonSimNhits[3];
   Int_t           fProtonSimCharge[3];
   Int_t           fProtonSimID;
   Int_t           fProtonSimPDG;
   Int_t           fProtonSimType;
 TVector3        *fDeltaRecPt;
   Double_t        fDeltaRecPtfX;
   Double_t        fDeltaRecPtfY;
   Double_t        fDeltaRecPtfZ;
 TVector3        *fDeltaSimPt;
   Double_t        fDeltaSimPtfX;
   Double_t        fDeltaSimPtfY;
   Double_t        fDeltaSimPtfZ;
   Double_t        fDeltaRecPhi_t;
   Double_t        fDeltaSimPhi_t;
 TVector3        *fRecKmu_t;
   Double_t        fRecKmu_tfX;
   Double_t        fRecKmu_tfY;
   Double_t        fRecKmu_tfZ;
 TVector3        *fSimKmu_t;
   Double_t        fSimKmu_tfX;
   Double_t        fSimKmu_tfY;
   Double_t        fSimKmu_tfZ;
 TVector3        *fRecKp_t;
   Double_t        fRecKp_tfX;
   Double_t        fRecKp_tfY;
   Double_t        fRecKp_tfZ;
 TVector3        *fSimKp_t;
   Double_t        fSimKp_tfX;
   Double_t        fSimKp_tfY;
   Double_t        fSimKp_tfZ;
   Double_t        fReconRecNuEnergy;
   Double_t        fReconSimNuEnergy;
   Double_t        fReconRecTgtMass;
   Double_t        fReconSimTgtMass;
 TLorentzVector  *fNeutronToy;
   TVector3        fNeutronToyfP;
   Double_t        fNeutronToyfE;
 TVector3        *fNeutronToyPt;
   Double_t        fNeutronToyPtfX;
   Double_t        fNeutronToyPtfY;
   Double_t        fNeutronToyPtfZ;
   Double_t        fNeutronToyPtMag;
   Double_t        fNeutronToyAlpha;
   Double_t        fNeutronToyQ2;
 TVector3        *fNeutronToyVertex;
   Double_t        fNeutronToyVertexfX;
   Double_t        fNeutronToyVertexfY;
   Double_t        fNeutronToyVertexfZ;
 TVector3        *fNeutronToyEndPos;
   Double_t        fNeutronToyEndPosfX;
   Double_t        fNeutronToyEndPosfY;
   Double_t        fNeutronToyEndPosfZ;
 TVector3        *fNeutronToyFlightPath;
   Double_t        fNeutronToyFlightPathfX;
   Double_t        fNeutronToyFlightPathfY;
   Double_t        fNeutronToyFlightPathfZ;
   Double_t        fNeutronToyChi2;
   Int_t           fNeutronToyNDOF;
   Int_t           fNeutronToyNdets[2];
   Int_t           fNeutronToyNhits[3];
   Int_t           fNeutronToyCharge[3];
   Int_t           fNeutronToyID;
   Int_t           fNeutronToyPDG;
   Int_t           fNeutronToyType;
 TLorentzVector  *fProtonToy;
   TVector3        fProtonToyfP;
   Double_t        fProtonToyfE;
 TVector3        *fProtonToyPt;
   Double_t        fProtonToyPtfX;
   Double_t        fProtonToyPtfY;
   Double_t        fProtonToyPtfZ;
   Double_t        fProtonToyPtMag;
   Double_t        fProtonToyAlpha;
   Double_t        fProtonToyQ2;
 TVector3        *fProtonToyVertex;
   Double_t        fProtonToyVertexfX;
   Double_t        fProtonToyVertexfY;
   Double_t        fProtonToyVertexfZ;
 TVector3        *fProtonToyEndPos;
   Double_t        fProtonToyEndPosfX;
   Double_t        fProtonToyEndPosfY;
   Double_t        fProtonToyEndPosfZ;
 TVector3        *fProtonToyFlightPath;
   Double_t        fProtonToyFlightPathfX;
   Double_t        fProtonToyFlightPathfY;
   Double_t        fProtonToyFlightPathfZ;
   Double_t        fProtonToyChi2;
   Int_t           fProtonToyNDOF;
   Int_t           fProtonToyNdets[2];
   Int_t           fProtonToyNhits[3];
   Int_t           fProtonToyCharge[3];
   Int_t           fProtonToyID;
   Int_t           fProtonToyPDG;
   Int_t           fProtonToyType;
 TVector3        *fDeltaToyPt;
   Double_t        fDeltaoToyPtfX;
   Double_t        fDeltaToyPtfY;
   Double_t        fDeltaToyPtfZ;
 TLorentzVector  *fProtonDet;
   TVector3        fProtonDetfP;
   Double_t        fProtonDetfE;
 TVector3        *fProtonDetPt;
   Double_t        fProtonDetPtfX;
   Double_t        fProtonDetPtfY;
   Double_t        fProtonDetPtfZ;
   Double_t        fProtonDetPtMag;
   Double_t        fProtonDetAlpha;
   Double_t        fProtonDetQ2;
 TVector3        *fProtonDetVertex;
   Double_t        fProtonDetVertexfX;
   Double_t        fProtonDetVertexfY;
   Double_t        fProtonDetVertexfZ;
 TVector3        *fProtonDetEndPos;
   Double_t        fProtonDetEndPosfX;
   Double_t        fProtonDetEndPosfY;
   Double_t        fProtonDetEndPosfZ;
 TVector3        *fProtonDetFlightPath;
   Double_t        fProtonDetFlightPathfX;
   Double_t        fProtonDetFlightPathfY;
   Double_t        fProtonDetFlightPathfZ;
   Double_t        fProtonDetChi2;
   Int_t           fProtonDetNDOF;
   Int_t           fProtonDetNdets[2];
   Int_t           fProtonDetNhits[3];
   Int_t           fProtonDetCharge[3];
   Int_t           fProtonDetID;
   Int_t           fProtonDetPDG;
   Int_t           fProtonDetType;
 TVector3        *fDeltaDetPt;
   Double_t        fDeltaDetPtfX;
   Double_t        fDeltaDetPtfY;
   Double_t        fDeltaDetPtfZ;
 TLorentzVector  *fNeutrinoSim;
   TVector3        fNeutrinoSimfP;
   Double_t        fNeutrinoSimfE;
 TVector3        *fNeutrinoRec;
   Double_t        fNeutrinoRecfX;
   Double_t        fNeutrinoRecfY;
   Double_t        fNeutrinoRecfZ;
 TVector3        *fNeutrinoGuess;
   Double_t        fNeutrinoGuessfX;
   Double_t        fNeutrinoGuessfY;
   Double_t        fNeutrinoGuessfZ;
   Double_t        fdNeutrinoGuessNeutrinoSim;
   Double_t        fdNeutrinoGuessMuonRec;
   Double_t        fGuessCorrection;
   Double_t        fMuonPtGuess;
   Double_t        fOsciMuMuProb;
   Int_t           fkOsciMuSurvive;
   Int_t           fNeutrinoParentPDG;
   Int_t           fNeutrinoParentType;
   Int_t           fNeutrinoType;
   Int_t           fNeutMode;
 TVector3        *fNeutrinoParentDecPointSim;
   Double_t        fNeutrinoParentDecPointSimfX;
   Double_t        fNeutrinoParentDecPointSimfY;
   Double_t        fNeutrinoParentDecPointSimfZ;
   Int_t           fMultiplicity;
   Int_t           fMuonNSim;
   Int_t           fMuonNRec;
   Int_t           fProtonNSim;
   Int_t           fProtonNRec;
   Double_t        fMeanBetaRec;
   Double_t        fMeanGammaRec;
   */